Water discharge variations control fluvial stratigraphic architecture in the Middle Eocene Escanilla formation, Spain

Ancient fluvial deposits typically display repetitive changes in their depositional architecture such as alternating intervals of coarse-grained highly amalgamated (HA), laterally-stacked, channel bodies, and finer-grained less amalgamated (LA), vertically-stacked, channels encased in floodplain deposits. Such patterns are usually ascribed to slower, respectively higher, rates of base level rise (accommodation). However, “upstream” factors such as water discharge and sediment flux also play a potential role in determining stratigraphic architecture, yet this possibility has never been tested despite the recent advances in the field of palaeohydraulic reconstructions from fluvial accumulations. Here, we chronicle riverbed gradient evolution within three Middle Eocene (~ 40 Ma) fluvial HA-LA sequences in the Escanilla Formation in the south-Pyrenean foreland basin. This work documents, for the first time in a fossil fluvial system, how the ancient riverbed systematically evolved from lower slopes in coarser-grained HA intervals, and higher slopes in finer-grained LA intervals, suggesting that bed slope changes were determined primarily by climate-controlled water discharge variations rather than base level changes as often hypothesized. This highlights the important connection between climate and landscape evolution and has fundamental implications for our ability to reconstruct ancient hydroclimates from the interpretation of fluvial sedimentary sequences.

. Conceptual figure explaining fluvial architecture as a function of the Accommodation (A) to Sediment supply (S) ratio and water discharge (Q w ). (a) Low A/S results in the deposition of high gradient fluvial channels with a high degree of amalgamation and an overall progradation of the system. Under high A/S, low gradient fluvial channels with a low degree of amalgamation are deposited with an overall retrogradation of the system. Figure based on the ideas by e.g., 4,10 (b) As Q w increases, low gradient channels with a high degree of amalgamation are deposited with an overall progradation of the system while as Q w decreases, high gradient channels with a low degree of amalgamation are deposited with an overall retrogradation of the system. In this scenario, base level rise is constant and changes in fluvial architecture are solely driven by changes in Q w .
Stratigraphic arrangement. This study focuses on 3 fining-upward sequences at Olson. Each finingupward sequence consists of a High Amalgamation (HA) package which can be defined as a thick (5-12 m) and laterally extensive (~ 600-2000 m) channel-body complex (90% channel to 10% floodplain) typically consisting of multiple stories (1)(2)(3)(4), and a Low Amalgamation (LA) interval which can be defined as a floodplain dominated (30% channel to 70% floodplain) interval with isolated channel bodies that are 2 -4 m thick and approximately 100 -500 m laterally extensive and typically consist of a single storey, although multiple stories can be found as well. On the field, the thickness of each such fining-upward sequence is around 35 -45 m and a HA interval typically lies above a thick (~ 5 -7 m) floodplain from the 'LA' interval of the underlying sequence (Fig. 3). Several different outcrop photographs can also be seen in Fig. 4 35 .

Results
Grain size, bankfull depth and palaeoslope evolution. Coarse grain size fraction (> 4 mm) (Fig. 5) was measured at 180 stations distributed across the studied stratigraphic interval. Coarse grain size fraction in channel bodies has an average range of grain sizes from 6 ± 2 [mm] to 30 ± 2 [mm] over the studied stratigraphic interval (Fig. 5). At the scale of individual sequences, HA intervals have a grain size of 19 ± 1 [mm] (average value ± standard error, N = 77), while LA intervals have a grain size of 17 ± 1 [mm] (N = 103), indicating that channel body grain size is 2 mm (~ 10%) larger in channels with a high degree of amalgamation. This relatively small difference is nevertheless statistically significant at the 95% confidence level (t value = 2.98, p value = 0.003, power = 0.91, dof = 178). It is noteworthy that the bulk grain size of LA intervals would be significantly lower if the fine-grained sediments of the floodplain were taken into consideration (50 -70% of LA intervals).
Bankfull depths based on preserved storey thickness (Fig. 5) reveal a trend of higher depths in the HA intervals and a substantial decrease in the LA intervals. HA intervals have a depth of 5.0 ± 0.4 [m] (average value ± standard error, N = 45), while LA intervals have a depth of 3.5 ± 0.3 [m] (N = 49), i.e., an increase of 40% of flow depth during deposition of the HA intervals. These field observations suggest that the palaeohydrology of channels comprising the HA and LA intervals is not the same. A t test on bankfull depth data rejects the null hypotheses that HA and LA intervals have the same average values at the 95% confidence level (t value = 6.83, p-value = 0.8e −9 , power = 0.67, dof = 92).
Palaeoslope estimates, obtained using the equation proposed by Trampush et al. 31 (Eq. 1), are consistently lower in the HA intervals, and markedly higher into the LA intervals (Fig. 5). Palaeoslope estimates based on averaged field grain size and channel depth data of 7 storeys within HA stratigraphic intervals (77 grain size sampling stations wherein 100 -200 clasts were counted per station, and 45 channel depth estimates) have a palaeoslope of 5 × 10 −4 ± 5 × 10 −5 [m/m] (average value ± standard error, N = 7), equivalent to 0.03°, while 11 storeys within LA stratigraphic intervals (103 grain size sampling stations and 49 channel depth estimates) have a palaeoslope of 8 × 10 −4 ± 6 × 10 −5 [m/m] (N = 11), equivalent to 0.05°, and representing a 60% increase in slope in the LA interval. A t test on palaeoslope estimates (t value = −4.02, p value = 0.001, power = 0.16, dof = 15) rejects the null hypotheses that HA and LA intervals have the same average values at the 95% confidence level. It is important to remark that the data in Fig. 5 show considerable variability typical of a natural fluvial system, but which nevertheless highlights systematic differences between HA and LA stratigraphic patterns. www.nature.com/scientificreports/ Channel-belt width estimates and channel body geometry. Channel-belt widths are necessary to estimating total discharge and sediment flux (e.g., 29,44 . Despite excellent outcrop conditions, fully preserved channel cross-sections required to measure width in the field, are nevertheless difficult to assess due to their size and the 3D nature of the outcrops. To circumvent this limitation, we estimate channel-belt widths Both estimates are smaller than the "geobody" widths estimated by Labourdette & Jones 34 , and we hence consider them to represent adequate conservative values. Lastly, a cross plot of estimated channel belt depth and width in Olson primarily implies sheet-like channel geometry during HA intervals, while channel geometry during LA intervals, even though corresponds to sheet-like geometry, is more closer to the ribbon fluvial style with decreasing depth and width ( Fig. 6).
Palaeohydrology and sediment transport. ). This represents a 1.5-fold increase in bedload sediment flux during HA intervals (Fig. 7). Considering the additional preservation of floodplain material during LA intervals, these results suggest a significant increase of sediment transport during the highly amalgamated HA intervals, and also predict commensurate export of this clastic material downstream of the Escanilla fluvial system (e.g. to shallow-marine environments). www.nature.com/scientificreports/  www.nature.com/scientificreports/

Discussion
These results provide new insights into how palaeoslopes and palaeohydrology of the Middle Eocene Escanilla system evolved with alluvial channel stratal architecture. We find that the channel-belt palaeoslope is systematically lower in HA intervals and higher in LA intervals, and that this corresponds to greater (resp. lower) water discharge and bedload sediment flux during HA (resp. LA) intervals (Fig. 7). Moreover, unit discharge estimates, which are width independent (i.e., independent of the estimated channel-belt widths), also document a systematic increase of streamflow during HA intervals, and decrease during LA intervals ( Figure S5). Currently, the hypothesis to explain the observed fluvial architecture in Olson is that it results from base level (accommodation variations) driven by local/regional tectonics in the Ainsa piggy-back basin 34 . However, such hypothesis (1) requires a, yet to be documented, physical mechanism that would account for tectonic oscillations with a periodicity of the order of 10 5 years to explain the observed stratigraphic changes (Fig. 2c), and (2) (3) it would also imply an opposite relationship between aggradation and slope (Fig. 2a), i.e. alluvial slope would decrease with base level rise during LA intervals (and conversely, river bed slope would increase during periods of low accommodation). Instead, our observations taken together suggest that the stratigraphic architecture in this alluvial depositional zone is rather controlled by upstream factors which modulate water discharge variations, sediment transport, alluvial style and stacking pattern (Fig. 7). We discuss below the possible origin of such variations. Upstream forcing of water and sediment fluxes can be associated to tectonic and/or climatic changes, which are the two allogenic controls primarily influencing depositional systems in the continental domain [e.g., [51][52][53] ]. Tectonic activity mainly influences sediment flux, normally without altering water discharge, except in case of structural disruption/modification of the drainage network in the catchment zone. Climatic changes, on the other hand, primarily affect precipitation, evaporation, overland and channel flow, which all influence the erosive and transport capacity within the drainage system, and thus determine the calibre and load of the fluvial sediment flux 53,54 . The effect of climatic changes on the production and release of sediment from source areas has been modelled in many studies [e.g., 52,[55][56][57], and a link between climatic oscillations and fluvial architecture through variations in water discharge, sediment supply, and sediment size has also long been proposed in numerous works [e.g., [58][59][60][61] . Thus, we postulate that a plausible explanation for the fining-upward trend of each HA-LA sequence could be a result of periodically decreasing precipitation (high during HA) and associated sediment discharge. With decreasing channel discharge, the fluvial equilibrium profile would be forced towards a greater gradient (as in Fig. 2b), driving fluvial aggradation, floodplain preservation and vertical stacking of channels (LA intervals).

Figure 7.
Palaeoslope, total water discharge and total bedload sediment flux estimates. Palaeoslope evolution along with estimates of total water discharge and total bedload sediment flux relative to each sampled storey within the HA and LA intervals of the studied sequences. Black vertical bars denote the average value in HA intervals while grey bars denote the average value in LA intervals. It is important to note the relationship and cyclical pattern, as shown by the overall average values in HA and LA intervals, between the three parameters such that river slope is lower when total discharge and total flux are higher (and conversely). Unit fluxes (water and sediment) have similar patterns of variations and are provided in supplementary figures S5 and S6. www.nature.com/scientificreports/ On the contrary, during periods of enhanced channel discharge, the river system would be forced toward a lower equilibrium bed slope (Fig. 2b), driving lateral mobility, amalgamation of channels and inhibiting vertical aggradation (HA intervals). At present, it remains unclear exactly what climatic conditions prevailed during the periods when we document larger (respectively lower) channel discharge. Indeed, the estimated larger discharges within HA intervals do not necessarily correspond to more humid conditions because drier climates with enhanced seasonality can also lead to channel-forming discharge of larger magnitudes. This has been suggested for similar situations during the PETM in Spain [e.g., Chen et al., 2018 62 ] and in the Piceance basin, Colorado, USA, where Barefoot et al. 63 documented increased channel mobility associated with enhanced seasonality. Further work, for instance using pedogenic and geochemical climatic proxies, will be needed to explore this issue. The suggestion that climate has an important influence on stratigraphic architecture in the Escanilla system is consistent with the work of Armitage et al. 6 on a numerical modelling study applied to the Escanilla sediment routing system. These authors predicted increased sediment flux due to increased precipitation in the catchment area during the Middle Eocene Climatic Optimum (MECO), approximately around 40 Ma near the C18r to C18n transition, which corresponds roughly to the "Olson sheet", i.e., the HA interval of sequence 2. We indeed document an increase in water discharge and sediment flux at that time, which could also correlate well with the pulse of deltaic progradation recently documented by Peris Cabré et al. 64 in coeval MECO successions of the adjacent Jaca basin to the West.

What drove climate variations?
The repetitive pattern of the three fining-upward sequences studied here points at a recurrence in the controlling forcings at their origin. Moreover, other studies in the area have described several similar sequences above and below those studied in the present work [33][34][35] , suggest that this cyclic control may have persisted over an extended period of more than 2.8 Ma according to current stratigraphic constraints. Cyclical oscillations of the orbital configuration of the Earth-Sun system over so-called Milankovitch periodicities constitute a plausible mechanism to account for the observed cyclic changes in water discharge in the Escanilla fluvial succession. Milankovitch cyclicity is known to have played a major role in driving climate change during Earth's history [e.g., 65,66 ], in the Cenozoic [e.g., 67 ], and has also been invoked in several studies in the Spanish Pyrenees as an important factor controlling depositional sequences [e.g., 68 ].
To test this hypothesis, we tentatively calibrated the eccentricity model of Laskar et al. 69 with the succession in Olson by matching the base (at 40 m on Fig. 8) and top (at 210 m) of Chron C18n (magnetostratigraphic data of 42 to the base of Chron C18n.2n and top of Chron C18n.1n on the Geomagnetic Polarity Time Scale (GPTS 2020) (Fig. 8) 43 . The magnetozone C18n.1r has also been tentatively placed at 80-90 m, based on the presence of a thin reverse period yet with low/ intermediate quality paleomagnetic samples in the original dataset of Vinyoles et al. 42 (Fig. 8). Although preliminary, this first-order scaling suggests a correlation between higher discharge HA intervals and phases of increasing eccentricity (400 ky signal, marked with a thick orange outline on Fig. 8), and conversely, lower discharge LA intervals correspond to phases of decreasing 400 ky eccentricity (Fig. 8). www.nature.com/scientificreports/ Such a correlation is consistent with the physical rationale according which periods of eccentricity maxima promote larger precession/insolation amplitudes more prone to trigger high-magnitude, extreme, events 70 , while decreasing eccentricity would instead diminish climate variability (seasonality) 71 and the intensity of channelforming discharge events. This scenario is also coherent with the study of Kodama et al. 72 who demonstrated a systematic link between the 400 ky eccentricity maxima and increased terrigenous supply from fluvial source to coeval shallow-marine successions near Arguis in the neighbouring Jaca basin. Finally, such orbital control on stratigraphic cyclicity has also been shown to influence depositional style through control on siliciclastic sediment supply to the deep-marine Early-Middle Eocene Ainsa basin [73][74][75][76] and may thus have been an important factor throughout an extended period of time in the South-Pyrenean foreland despite the important tectonic activity linked with mountain building at that time.
Implications and conclusion. These findings have several fundamental implications for classical sequence stratigraphic predictions, sedimentary landscape evolution over geological timescales, the prediction of ancient hydroclimates during greenhouse forcing, and for industrial applications in resource exploration. First, our findings illustrate the dominant role of the sediment supply S on the stratigraphic architecture of a fluvial system, which contrasts with more classical hypotheses invoking a dominant role of accommodation A in determining fluvial stratal patterns [e.g., 7,10 ]. Second, our results suggest higher sediment fluxes during HA intervals, which is in agreement with classical sequence stratigraphic predictions in which such amalgamated packages are associated to low A/S conditions, and thus with greater down system transport of sediments to shallow-marine and deep-sea environments. Third, this work is also fundamental for our ability to reconstruct ancient hydroclimates from the sedimentary record and compare it to numerical model predicted response of the hydrologic cycle to warming conditions. For instance, modelling studies of the Palaeocene-Eocene Thermal Maximum (PETM) have suggested an intensification of the hydrological cycle on a global scale in response to greenhouse gas levels 77 . Results from this study are also consistent with the findings of Barefoot et al. 63 where higher seasonality during the PETM increased lateral channel mobility. During the Middle Eocene greenhouse conditions of the Escanilla succession, HA intervals therefore record periods of increased lateral channel mobility linked with channel-forming discharge events (floods) of higher magnitude than LA intervals which present lower lateral channel mobility (lower magnitude channel-forming discharge events). The framework of palaeoslope and palaeohydrological reconstruction across three fining-upward sequences in the Middle Eocene aged Escanilla Formation, for the first-time, documents lower river slope during higher water discharge and sedimentary flux with the deposition of river channels having a high degree of amalgamation, and higher slope during lower water discharge and sedimentary flux with the deposition of river channels having a low degree of amalgamation. The studied fining-upward sequences together with our finding of cyclic variations in water discharge and sediment flux represents a major paradigm shift by suggesting climate may have controlled the entire sedimentary landscape evolution during the Middle Eocene greenhouse conditions instead of eustatic variations, and with palaeohydraulic reconstructions to test different options without the need of an independent eustatic curve.

Methods
Data collection. Palaeohydrological field data were collected sequence-wise in a lateral spatial domain and from the respective HA and LA intervals. Data from channel fill deposits included channel basal gravels for grain size distribution and storey thicknesses as flow depth estimates. This data along with uncertainties associated with individual measurements were propagated through a quantitative framework to reconstruct hydrological parameters such as flow depths, palaeoslope, flow velocities, water discharge rates, bedload sediment flux. A Google Earth panel with the location of sampling stations has also been provided in the supplementary material.
Grain size measurements were collected using the Wolman sampling procedure 78 . The longest axis was measured as a proxy for the intermediate b axis on 100-200 grains per sampling station. The procedure was performed on photographs taken with a Canon EOS 2000D camera of 24.1 Mpixels resolution on an outcrop area of 1 × 1 m 2 . Grains were measured at the nodal intersection of a virtual grid such that a repeat count of grains is avoided. Measurements were made using ImageJ2 version 2.3.0/1.53f. The data obtained is normalized using the (psi) scale, a logarithmic scale with base two, to perform statistical analyses and obtain the 50th percentile (D 50 ) of the grain size distribution. Flow depth estimates are based on preserved storey thicknesses, channel-plug, and bar-scale clinoform heights measured using a laser range finder (TruPulse model 200) and following the procedure outlined in Mohrig et al. 79 and Kelly 48 . It is important to note that while preserved thicknesses are lower than the original flow depths, preserved thicknesses do not severely underestimate the original depth 32,80 . Quantitative paleohydrology. Palaeoslopes were estimated using the empirical equation proposed by Trampush et al. 31 (Eq. 1). We use this equation as our grain size measurements in a few instances are less than the 8 mm threshold required to use the Shields stress inversion approach 32 .
It is an empirical equation, motivated by theoretical considerations, and provides a relationship between the channel slope (S), median grain size ( D 50 ), and bankfull depth ( H bf ). α 0 , α 1 and α 2 are three empirical coefficients with values of − 2.08 ± 0.0015 (mean ± standard error), 0.2540 ± 0.0007 and − 1.0900 ± 0.0019, respectively. An average palaeoslope (± standard error) value has been estimated, per interval, using average median grain size (1) log S = α 0 + α 1 log D 50  www.nature.com/scientificreports/ values and average bankfull depths along with their respective standard errors. Average palaeoslope estimates are presented in [m/m], for example, a palaeoslope value of 0.001 represents aggradation of 1 m per 1000 m. Channel-belt width W can be estimated using empirical scaling relations when direct measurements are not possible on the field. We estimate channel-belt widths using the relationship, W = 8.8H bf 1.8245 . Where possible, channel plug widths were estimated in the field using a laser range finder (TruPulse model 200) while widths from lateral accretion deposits was estimated using the procedure outlined in Greenberg et al. 44 .
Flow velocity, U was calculated using Manning's equation (Eq. 2) where n = 0.03 ± 0.005 is the Manning's coefficient, R is the hydraulic radius approximated by channel flow depths and S is slope. Flow velocity data can be found in the supplementary material.
Total water discharge Q w was calculated using Eq. (3) For unit water discharge, W = 1. Total bedload sediment flux Q s was calculated using the Meyer-Peter and Muller equation (Eq. 4).
Effectiveness of palaeohydrological reconstructions. While palaeohydrological information can be extracted for ancient systems, there remain concerns over its accuracy 29 . To address these concerns, we evaluate a few limitations in this study. Firstly, there is a possibility of under-estimating true depths based on storey thickness due to the incomplete preservation and erosion by the overlying storey. An under estimation of flow depth would consequently affect the slope, flow velocity, discharge, and flux estimates. Secondly, our total discharge and flux estimates are based on channel-belt widths. While individual channel widths can easily be determined for modern rivers, the inherent problem in doing the same for ancient systems, apart from low preservation, lies in determining the number of active channels. This could significantly affect our estimates such that a high number of active channels would imply significantly higher discharge and sediment transport rates.
Statistical tests. Uncertainty on results reported in this study consist of the standard error of the mean (SE) calculated as = SD √ n , where SD is the standard deviation and n is sample number. Uncertainty propagation was carried out using the uncertainties package on Python (Spyder version 4.0.1). Statistical analyses were performed on Python (Spyder version 4.0.1). To check for data normality, the Shapiro-Wilk test was performed using the 'scipy.stats.shapiro' package. The null hypothesis that the data is normally distributed cannot be rejected when the p value is greater than 0.05 at the 95% confidence level. To check for statistical significance, a two-sided t test was performed for normally distributed data using the 'scipy.stats.ttest_ind' package for the null hypothesis that two independent samples have an identical average value. For non-normally distributed data, a Kruskal-Wallis test was performed using the 'scs.kruskal' package for the null hypothesis that the median value of all groups is similar. The null hypothesis can be rejected when the p-value is less than 0.05 at the 95% confidence level. The degree of freedom (dof) was estimated as '(nx + ny) − 2' where 'nx' and 'ny' are the lengths of the two independent parameters. Power analysis of t tests was performed, using 'pingouin.power_ttest2n', to detect Type II errors.
Pingouin is an open-source package written in Python 3 81 .

Data availability
All data generated and analysed in this study, including locations of sampling stations have been provided as supplementary material files accompanying this manuscript.